function [output] = SK_main(input_array)

% data input
method = input_array.method;
data = input_array.data;
Xcoord = input_array.Xcoord;
Ycoord = input_array.Ycoord;
alt_vec = input_array.alt_vec;

D = data.D;
X = data.X;
Y = data.Y;
X_mat = [D X];
        
switch lower(method)
    case 'ols'
        betaHat_ols = X_mat\Y;
        beta0_hat = betaHat_ols(1); % OLS estimate
    case 'iv'
        % IV method
        Z = data.Z;
        Z_mat = [Z X];
        beta_iv = (Z_mat'*X_mat)\(Z_mat'*Y);
        beta0_hat = beta_iv(1);
end

p_null = test_harSunKim(method,data,Xcoord,Ycoord);

p_alt = 0*alt_vec;
for i_alt = 1 : length(alt_vec)
    data_alt = data;
    data_alt.Y = data.Y+data.D*alt_vec(i_alt);
    p_alt(i_alt) = test_harSunKim(method,data_alt,Xcoord,Ycoord);
end


%% output
output.estimate = beta0_hat;
output.p_null = p_null;
output.p_alt = p_alt;

